Satisfying the Fluctuation Theorem in Free Energy Calculations with Hamiltonian 

Replica Exchange 
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A novel error measure, referred to as the hysteresis error, is developed from the Crooks fluctuation 
theorem to evaluate sampling quality in free energy calculations. Theory and numerical free energy 
of hydration calculations are used to show that Hamiltonian replica exchange provides a direct 
route for minimizing the hysteresis error. Replica exchange swap probabilies yield the rate at which 
the hysteresis error falls with simulation length, and this result can be used to decrease bias and 
statistical errors associated with free energy calculations based on multicanonical simulations. 

PACS numbers: 



I. INTRODUCTION AND OVERVIEW 

Free energies of solvation provide quantitative assess- 
ments of driving forces for spontaneous processes such 
as protein folding, binding, self-assembly, and solubility. 
Formally, the free energy of solvation in the canonical en- 
semble is the free energy change AF associated with the 
transfer of a solute from the gas phase to a fixed position 
in the solvent [ij . Operationally, one has access to a range 
of techniques to obtain estimates for AF Kirkwood 
showed that one could introduce arbitrary parame- 
ters into potential functions and continuously vary the 
degree of coupling between specific molecules in a dense 
fluid. The device of coupling parameters leads to simple 
expressions for chemical potentials of any component of 
the fluid. If the component is the solute molecule, which 
is transferred from the gas phase into the solvent, then 
a single coupling parameter A, where < A < 1, modu- 
lates solute-solvent interactions in the system's potential 
function. The limits A = and A = 1 correspond to the 
pure solvent and solvent plus fully grown solute, respec- 
tively. Intermediate values of A correspond to potential 
functions that include only a part of the solute-solvent 
interactions. The Kirkwood coupling parameter plays a 
central role in equilibrium methods for calculating A_F. 
One carries out a series of independent canonical simula- 
tions where each simulation is associated with a distinct 
potential function, characterized by a specific A value. 
As it samples the equilibrium ensemble, each simulation 
generates a series of work values, which are then used 
to estimate the free energy change across the entire A 
schedule. 

The multicanonical approach described above takes 
advantage of the simple formalism developed by Kirk- 
wood for calculating AF. However, in practice, standard 
free energy calculations based on multicanonical simula- 
tions are plagued by slow convergence and inaccurate es- 
timates of AF Q. Errors may be divided into statistical 
and bias (or finite sampling) errors Q . The former stem 
from the fluctuations of the free energy estimator, and 
can be estimated by block averaging or bootstrap meth- 
ods 0, Q. Since the statistical error decreases as the 



inverse square root of simulation length, it is frequently 
used as an indicator of the convergence of the multi- 
canonical simulation. While statistical errors are random 
fluctuations of short simulation results about some mean 
value, the bias error is an error of the mean value it- 
self, and it changes with simulation length. As discussed 
by Zuckerman and WoolfjoJ, bias errors have two causes: 
the free energy estimates are nonlinear averages; and, 
the work distributions on which such estimates are based 
will typically have long tails which are rarely sampled, 
and yet these are important to the average. The latter 
point is important: rare events dominate free energy es- 
timates, and one seldom observes these events in short 
simulations. As a result, the average drifts with simu- 
lation length, resulting in inaccurate estimates for AF 
from bias error even when the statistical error is small. 
The magnitude of the bias error is difficult to quantify 
directly, as it requires knowledge of the actual free en- 
ergy difference, the very quantity we wish to determine. 
Furthermore, small fluctuations in the estimate for AF 
may not be indicative of convergence, but rather of inade- 
quate sampling of the rare but important conflgurations. 
To address these problems, we develop an alternate mea- 
sure of free energy error, one based on deviations from 
equilibrium distributions. 

Crooks 10] derived a fluctuation theorem (appendix 
lA II) valid for stochastic, microscopically reversible dy- 
namics, which relates the distribution of dissipated work 
values along a forward and reverse path as. 
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Here, (3 = (ksT)-^, Pf[!3Wd) is the probability dis- 
tribution for dissipated work associated with switching 
A from Ao to Ai, and Pr{—(3Wd) is the corresponding 
distribution for the reverse process. If the canonical sim- 
ulations for each value of A sample the equilibrium en- 
semble adequately, then the distributions of dissipated 
work obtained over the course of free energy calculations 
will satisfy Eq. ([T]). 

In this work, we develop a readily measured error 
estimate, the hysteresis error eH^ which quantifles the 
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degree to which observed work distributions obey the 
Crooks fluctuation theorem. Hamihonian rephca ex- 
change, a multicanonical equilibration technique, effec- 
tively reduces the hysteresis error. We relate the average 
replica exchange swap probability to the degree of over- 
lap between equilibrium ensembles, as well as to the rate 
at which en falls. Based on this, we may construct an 
optimized A schedule to further minimize the hysteresis 
error for an entire simulation. 

The remainder of this presentation is organized as fol- 
lows: the theory section introduces the hysteresis error 
in the context of the Crooks fluctuation theorem followed 
by a formal illustration of how Hamiltonian replica ex- 
change minimizes ch] the definition of swap probability 
as a measure of the overlap between different equilibrium 
ensembles; and a connection between the amount of over- 
lap and minimization of e}j. We calculate the free energy 
of hydration for acetamide to demonstrate how to esti- 
mate en and minimize this error using replica exchange 
coupled to standard multicanonical simulations. We con- 
clude with a summary and a discussion of the features of 
our methodology. 

II. THEORY 
A. Background 

The free energy of replica i in the canonical ensemble 
at temperature T, whose potential t/i(r) = U{T,\i) is a 
function of system configuration T and the parameter Ai, 
is formally given as [ll[, 

= /3-iln|y"drexp[-/3i7,(r)]|. (2) 

At equilibrium, the probability of observing configuration 
r is given as, 

p,(r)=exp{/3[i^,-C/,(r)]}. (3) 

To calculate the free energy change 5F associated with 
switching the Hamiltonian from Uq to Ui we perform 
simulations at Ag and Ai, and calculate the forward and 
reverse work as, 

W^{T) = Ui{T) - C/o(r), (4a) 

w^^(r) = c/o(r) - c/i(r). (4b) 

For the forward and reverse work values the configura- 
tion r is typically drawn from the equilibrium ensemble 
of Uo and Ui, respectively. The Free Energy Perturba- 
tion (FEP) method [l^l utilizes forward and reverse work 
distributions to provide two independent estimators for 
SF, 

SF^EP = -/3-Mn(exp(-/3iy^)}o, (5a) 
SFPep = +/3-Mn(exp(-/3M^^))i, (5b) 



where the forward estimator SFp^p utilizes forward work 
values from the simulation at Uq, and the reverse estima- 
tor the reverse work from Ui. Note that in both cases 
SF is associated with the process of switching Aq — > Ai. 
These two estimators have different convergence rates . 
Therefore, while in practice the two estimates should be 
equal, in simulations with finite sampling they are gen- 
erally different. 

Another free energy estimator, the Bennett Accep- 
tance Ratio [31, uses both the and distribu- 
tions to obtain a free energy estimate. It is generally 
more accurate (l^ and is employed later in this paper 
for numerical free energy estimates, but will not be con- 
sidered for theoretical development. 



B. The Hysteresis Error 

The hysteresis error en is defined as the difference be- 
tween the forward and reverse SFfep estimates, 

eH = 5FfEp~SF^Ep. (6) 

en has contributions from both the statistical and bias 
error of the FEP estimators [1, Q. The bias error of 
the two estimators is typically in the opposite direction. 
While the statistical error may dominate the en for a 
given simulation, in averages over multiple short simula- 
tions the dominant contribution to the average hysteresis 
error is the sum of the forward and reverse FEP bias. 

We take en as a measure of sampling quality and aim 
to minimize its magnitude between all pairs of neigh- 
boring replicas. The validity of using en as a general 
sampling error is based on a relationship between it and 
the fluctuation theorem of Crooks(IT|), derived below. 

Switching the parameter Aq — > Ai (and vice versa) is 
equivalent to performing non-equilibrium work; the dif- 
ference between the work performed and the free energy 
change of the system is the dissipated work, defined in 
the forward and reverse direction as, 

wE{r)^w''{r)~SF, (7a) 

Wg{r) = T4^^(r) + SF. (7b) 

Crooks [Toj equates and Wj^ to the entropy produc- 
tion caused by changing Aq — > Ai and Ai — > Aq, respec- 
tively, for the given configuration. 

The distributions Pf{Wd) and Pr{Wd) give the prob- 
ability of realizing a specific value for the dissipated work 
in the forward and reverse directions, respectively. The 
distributions are related to each other by the fiuctuation 
theorem shown in Eq. ([T]), which we have re-derived in 
appendix IA II for the specific case of instantaneous switch- 
ing between configurations with different A values. In 
practice, Eq. ([T|) will not be satisfied exactly because of 
errors due to finite sampling. To take simulation errors 
into account, we rewrite Eq. ([T]) with an arbitrary error 
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term Cpj, and with observed (rather than ideal) dissi- 
pated work distributions Pp and P^ , 



PU^PWd)' 



(8) 



Eq. ^ is constructed such that the Crooks fluctuation 
theorem is recovered and e*prp = when the observed 
work distributions match the correct distributions. The 
hysteresis error en and the fluctuation error c^rp are re- 
lated to each other as, (see appendix lA 2[) . 



tH^-p ^ ln(exp(-/3e>5.))5, 



(9) 



where (•)* is defined as the average obtained from a fi- 
nite simulation. The more closely a simulation obeys the 
relationship ([T]), the smaller the hysteresis error eu, and 
vice versa. In the next section, we will discuss methods 
to reduce en, which in turn leads to the satisfaction of 
the Crooks fluctuation theorem. 



C. Replica Exchange 

In a Hamiltonian replica exchange (isl . [l6j simulation, 
Monte Carlo moves are employed to exchange configu- 
rations r (or equivalently, parameters A) between two 
replicas with the probability, 
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swap J 1 5 



where, 



AU. 



swap 



c/o(ri) + c/i(ro) 
-c/o(ro)-c/i(ri). 
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(lib) 
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Fq and Fi denote configurations drawn at random from 
the equilibrium ensembles of Uq and Ui , respectively. For 
convenience, we write 7 = (Fg, Fi) as a pair of such con- 
figurations, and 7' = (Fi,Fo) is the swapped configura- 
tion pair. 

Since Fq and Fi are independent configurations, we can 
consider the probability of sampling Fq in the equilibrium 
ensemble of Uq and sampling Fi in the equilibrium en- 
semble of C/i; this is the native probability Pn{i)- Anal- 
ogously, the joint probability of sampling the swapped 
configurations, Fi from po and Fq from pi is given as 
pW(7): 

PAr(7) = po(Fo)pi(Fi), (12a) 
p'^(7) = po(Fi)pi(Fo)=pAr(7')- (12b) 

Replica exchange swaps are conveniently visualized by 
plotting the independent configurations Fq and Fi along 
orthogonal axes and the equilibrium ensemble of the sys- 
tem as an isocontour of pisi, illustrated in Fig. [Ija). 



At equilibrium, the relative probability of observing a 
pair of replicas in their swapped versus native configura- 
tions is, 



Pn 



exp(-/3AC/s„ap), 



(13) 



which is derived with definitions (I12p , ([3]) and (|llap . We 
will refer to this as an inter-replica equilibrium relation- 
ship. 

In an infinitely long simulation, (|13p will be satisfied 
exactly, but this will generally not be the case for fi- 
nite simulations, where inadequate sampling of config- 
uration space will result in inaccurate probability esti- 
mates. However, in simulations with replica exchange 
we expect the inter-replica equilibrium relationship to be 
satisfied more closely than in simulations without replica 
exchange, because the swap move distributes configura- 
tion pairs in such a way as to satisfy Eq. (|13p . To il- 
lustrate, consider the system in Fig. [T]^a) where the Uq 
replica is presumed to be stuck in the left lobe of the po 
distribution because of a kinetic barrier. Without replica 
exchange, only the shaded region of pn will be sampled 
accurately. The simulation will not have a correct es- 
timate for p'j^ji'jc) — Pn{i'c)-, since po for the swapped 
configuration, never having been observed, will be in- 
accurate. Consequently, Eq. (|13p will not hold. Replica 
exchange directly populates swapped configurations (e.g., 
7^), thereby improving the statistics of p^ and allowing 
inter-replica equilibrium to be achieved more quickly for 
all configurations in piq. 

The degree to which Eq. (|13p is satisfied determines 
the magnitude of the hysteresis error. To illustrate this, 
suppose that the distribution p^ has some small error 
Pe(Fo, Fi) due to finite sampling, so that we write (p^ -I- 
Pe) as the numerator in Eq. (|13p . In appendix lA 31 we 
show, by integrating over all configuration pairs, that the 
relationship between the hysteresis error and the error of 
sampling the swapped distribution, p^ is, 



-r 



dForfFiPe 



(14) 



The hysteresis error, then, will be minimized when the es- 
timated swapped configuration probabilities p^ are con- 
sistent with the equilibrium distribution. Since replica 
exchange populates the swapped configurations directly, 
it provides an efficient route to minimizing eu- 



D. Swap Probability 

Analysis of the average swap probability is complicated 
by the fact that the Metropohs function (Eq. HU])) is not 
analytical. For the purposes of interpreting this quantity, 
we will instead consider the Fermi swap probability. 
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FIG. 1: A graphical representation of replica exchange, (a) The independent (high dimensional) configuration spaces To and 
Fi have probability distributions po and pi, respectively, and the joint equilibrium ensemble pjv is drawn over this domain. The 
po system has a kinetic barrier (represented by the two disconnected lobes) and with no replica exchange the system explores 
only the configurations of the shaded domain. A replica exchange swap is a refiection of the configuration pair 7 about the 
Fo = Fi diagonal axis, and three swap attempts are shown: the configuration pair 7a swaps successfully and becomes 7^, but 
it does not sample otherwise inaccessible regions; a swap of 7;, fails because 7^ is not in the equilibrium ensemble; and the swap 
of 7c succeeds and allows the system to explore otherwise inaccessible regions of phase space, (b) The equilibrium domain pjv 
and its swapped image p^ are drawn. Swaps are feasible only for configuration pairs which belong to both pjv and p^. This 
overlap region, labeled Pswap, is the domain where the integrand of Eq. (|17b[l is large, and its size corresponds to the average 
swap probability, (c) The overlap of the po and pi distributions along the common configuration Fo — Fi. For the hysteresis 
error to converge, the Ao simulation must observe configurations where pi > po, and the Ai simulation must adequately sample 
the region po > pi. The frequency with which this occurs is given by (psmap)- 



where f{x) is defined as, 

/(x) = l/[l + exp(x)]. (15) 

(See [il] for discussion). We use Pswap to denote the 
Fermi swap probability and Pswap for the Metropohs 
swap probabihty; while the theoretical development uses 
Pswap, replica exchange moves are accepted/rejected us- 
ing Pswap- A simulation with either the Metropolis or 
Fermi swap probability will yield a Boltzmann distri- 
bution of swapped and unswapped configurations (Eq. 
([T^ ) . While the exact numerical values of the Fermi and 
Metropolis swap probabilities will differ somewhat, their 
qualitative behavior and the conclusions drawn here will 
hold for both. 

The average Fermi swap probability for two systems 
evolving independently is, 

(Pswap) = {{f{f3AU swap)) 0)1, (16a) 

= J dTodVipNfiPAUswap), (16b) 

which can be written as, 

iPswap) = ((^^) ) , (17a) 

\\PN + Pn / a/ I 

= [dTodV,^^. (17b) 
J Pn + Pn 



The integrand of (jl7bp is a normalized probability of ob- 
serving a given configuration pair, and the average swap 
probability is then the overlap oi pN and p'jq. See Fig. 
[Tfb) for a graphical interpretation. Thus, a large aver- 
age swap probability implies a large overlap between the 
equilibrium distributions of the two replicas, and a low 
(Pswap) indicates that the configurations these replicas 
adopt are distinct. 

We can expand (|16ap in a Taylor series about A = 
+ S\. To leading order in Sx, we find that in the 
neighborhood of Ao the average swap probability is, (see 
appendix lA 4p , 

iPswap) ^ 2 - (18) 

where 

Ca, then, determines the rate at which the average swap 
probability declines as the difference in A between the 
two replicas, S\, increases, although this linear analysis 
is accurate only for small 5\. 
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E. Swap Probability and the Hysteresis Error 
Convergence Rate 

We now demonstrate that the average swap probabihty 
between two rephcas gives a measure of how quickly the 
hysteresis error decreases, on average, over the course 
of a simulation. The hysteresis error is the difference 
between the forward and reverse SFfep, and since the 
forward and reverse FEP estimators do not converge at 
equal rates @ , it is the slower of these which governs the 
convergence of en- 

We may rewrite Eq. ([5a| as, 



(19) 



For this to hold, we must sample configurations where 
< 0; since the dissipated work is on average greater 
than zero by the second law of thermodynamics, such 
configurations tend to be rare [13. As a result, the 
convergence rate of SFp^p is governed by the probabil- 
ity of observing negative dissipated forward work values. 
Likewise, the convergence of SFp^p is dictated by ob- 
servations of < 0. We can understand this criterion 
graphically with the relationships, (see appendix I A ip . 



Po(ro) 

Pi(ro) 
/>i(ri) 
Po(ri) 



(20a) 
(20b) 



In the context of Fig. [He), observing < corre- 
sponds to sampling configurations from the po distribu- 
tion where pi > po, and for W§ < we require po > pi 
when sampled from the pi distribution. 

Turning our attention to the average swap probabil- 
ity, we note that /^Uswap, which is the sum of and 
Wj^, is negative whenever p'^ > pjv (by Eq. (fT^ ). Con- 
figurations for which this is the case are sampled by a 
simulation only in the lower-right half of the domain la- 
beled Pswap in Fig- HJb). The larger this domain, whose 
size is given by the average swap probability, the more 
frequently negative values of and W§ are observed, 
and the more quickly the hysteresis error converges. A 
numerical confirmation of this argument, that low swap 
probabilities correspond to large hysteresis errors and 
vice versa, is demonstrated in the results section. 



III. METHODS 

The computational system consists of 21 replicas, each 
with a different A, which are simulated independently 
to obtain equilibrium statistics. The parameter A con- 
trols the non-bonded interactions between an acetamide 
(ACE) solute and the water molecules. Two indepen- 
dent sets of simulations were performed, with and with- 
out replica exchange, in order to investigate the effect of 
this technique. 



The Lennard- Jones and Coulomb interactions between 
the water and ACE molecules are scaled by Xlj and Ac, 
respectively. We scaled both parameters simultaneously, 
such that Aij = Ac; the single parameter A then refers 
to both terms. The specific way in which the Lennard- 
Jones and Coulomb terms scale with A is described in 
appendix [Bl A varies across the 21 replicas from to 1 
in increments of 0.05. 

Each replica consists of 343 water molecules and one 
ACE molecule, which is rigid and whose position is fixed 
in the central box. All simulations were performed at 
constant temperature (298K) and volume (21. SA cubic 
box) using Metropolis Monte Carlo sampling. Parame- 
ters from the OPLS-AA force field [3 and 4-site TIP4P 
water model [3 were used to model the solute and sol- 
vent, respectively. Minimum image boundary conditions 
and spherical cutoffs were employed for the Coulomb and 
Lennard- Jones potentials. The cutoff radius was 10. 5A 
for electrostatic interactions and lOA for van der Waals 
interactions. Cutoffs were group-based for the former, 
and atom-based for the latter. No long-range correc- 
tions were employed. All simulations were carried out 
using the MCCCS Towhee [13] Monte Carlo simulation 
package [32j. 

The initial configurations for all replicas were identical 
and correspond to the end-point of a pre-equilibration 
run with ACE in water. For each replica, simulations 
consisted of 2 million cycles, where a cycle corresponds 
to 343 Monte Carlo moves; each move combines rotations 
and translations of a randomly chosen individual water 
molecule. The initial 10^ cycles were discarded for equili- 
bration. The average acceptance rate for all replicas was 
31%. 

The replica exchange simulation consists of a number 
of simulation rounds, where each replica evolves inde- 
pendently, separated by swap rounds, when a number 
of swap attempts take place. The length of the simula- 
tion round was drawn from a normal distribution with 
a mean of 500 and standard deviation of 50 cycles. 500 
cycles is the approximate energy autocorrelation "time" . 
The swap round consists of 21^ swap attempts between 
randomly selected replica pairs. Allowing swaps beyond 
neighboring replicas increases the efficiency of replica ex- 
change, by allowing a replica to traverse the entire range 
of A from to 1 more quickly than if only neighbor swaps 
were permitted [21]. 

During the course of the simulation, the native {Ui{Ti)) 
and foreign {Uj^i{Ti)) potential energies, as well as 
values for dU/dXc and dU/dX^j (where dU/dX — 
dU/dXLj + dU/dXc), were saved every 10 cycles. These 
were then post-processed to obtain the free energies, the 
hysteresis error, swap probabilities, and Ca, regardless of 
whether actual replica exchange swaps took place. The 
total free energy of hydration, AF, is the sum of all free 
energy changes {SF)i between neighboring replicas i and 
i + 1, calculated using the Bennett Acceptance Ratio 
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method [13| . 



M-l 



be reduced by decreasing the A-distance between them, 
and an optimal A schedule can reduce it for an entire 
simulation. 



where M is the total number of replicas. Similarly, the 
RMS hysteresis error crms is the root-mean-square of 
the hysteresis error (e//)i between neighboring replicas, 



f-RMS 



\ 



M~l 



Statistical errors for A_F were estimated using the 
bootstrap method With the simulation dataset con- 
sisting of N observations, we drew n* observations at 
random and with replacement to create one bootstrap 
estimate, AF*. This process was repeated 10,000 times, 
and the standard deviation among all the AF* is the 
estimated error of A_F. n* is the expected number of in- 
dependent observations in the dataset; here, n* = 1900 
with the assumption that there is one independent obser- 
vation per two internal energy autocorrelation "times" 



B. Hysteresis Error and Replica Exchange 

For a fixed A schedule, the hysteresis error may be 
reduced with either an improved sampling methodology 
like replica exchange, or longer simulations per replica. 
The effects of both approaches are illustrated in Fig. [2l 

Panel (a) shows eu for each neighboring replica pair. 
The hysteresis error is not uniform across all pairs, with 
spikes in the region A = 0.1 — 0.3. Replica exchange 
systematically reduces the hysteresis error for all pairs of 
replicas. 

Panel (b) illustrates how both longer sampling and 
replica exchange affect the hysteresis error. Block averag- 
ing shows that the average RMS hysteresis error declines 
consistently with longer simulations. This reduction can 
be improved with replica exchange; in fact, a simulation 
with replica exchange will achieve the same magnitude 
of RMS hysteresis error 5 times more quickly than one 
without replica exchange. 



IV. RESULTS 
A. Acetamide Free Energy of Hydration 

The hydration free energies we calculate for acetamide 
are in line with results obtained by other researchers, 
as shown in Table IH All numerical results differ some- 
what from experimental values due to differences in force 
field parameters. Our calculations were carried out in the 
canonical ensemble. Therefore, we obtain estimates for 
the Helmholtz free energy Af , whereas the experimental 
and other computational values obtain estimates for the 
Gibbs free energy, AG. However, the distinction between 
these two values should be negligible [1^. The consis- 
tency between our results and those of others serves to 
verify our implementation and sampling technique. 

TableUshows differences between results obtained with 
and without replica exchange. As expected from our the- 
oretical considerations, we find that the RMS hysteresis 
error is lowered by an order of magnitude when replica 
exchange is coupled to the multicanonical sampling pro- 
tocol. However, it should be noted that the statistical er- 
ror estimated using bootstrap remains unaffected. This 
is not an artifact of the bootstrap method used to esti- 
mate statistical errors. Instead, fiuctuations in estimates 
for 6F originate in fluctuations of the underlying work 
distribution, shown in Eq. ([1]). So long as both simula- 
tions sample the work distribution adequately, they will 
have similar statistical error associated with them. As a 
cautionary note, low statistical errors can also be caused 
by inadequate sampling of the appropriate work distri- 
butions. The statistical error between two replicas can 



C. Average Svirap Probability 

Fig. [3] shows downward spikes in the swap probabil- 
ity for values of A where the hysteresis error is large in 
Fig. mja). These results are consistent with the proposal 
that swap probability between two replicas is an indi- 
cator of the rate at which €h is minimized. The same 
region is characterized by a positive spike in C\, which 
is expected based on the relationship between the swap 
probability and C\ in Eq. (|18p . However, while the swap 
probability calculation requires two separate simulations, 
estimates of C\ can be obtained from just one. More- 
over, {pswap) varies as the distance between the replicas 
changes, complicating the interpretation if the A schedule 
is not uniform. Evaluation of C\ as a function of A us- 
ing a preliminary, coarse A schedule can identify regions 
where the swap probability is expected to be low, and can 
be used to construct optimal A schedules, as discussed in 
Sec. WM 

V. DISCUSSION 

A. Physical Interpretation of C\ Profile 

To gain a physical interpretation of the profile for Cx 
shown in Fig. [31 we plot in Fig. [l] the average water 
density in a 2.5A sphere surrounding the carbonyl car- 
bon of acetamide. The plot shows that water occupancy 
around the growing solute decreases rapidly in the range 
of A ~ 0.15. The expulsion and rearrangement of water 
molecules during cavitation leads to a large shift in the 
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(a) Acetamide Free Energy of Hydration: Current Work 





AF (kcal/mol) 


eRMS (kcal/mol) 


No Replica Exchange 
Replica Exchange 


-8.35 ± 0.051 
-8.14 ± 0.053 


0.120 
0.023 


(b) Acetamide Free Energy of Hydration: Literature 




AG (kcal/mol) 


Details 


MacCallum and Tieleman [24] 
Shirts et al. [25] 
Chang et al. [26] 
Udier-Blagovic et al. [27] 


-8.25 ± 0.26 
-8.20 " ± 0.03 
-8.54 ± 0.1 - 0.3 
-9.65 ± 0.3 - 0.5 


TIP4P, TI 
TIP3P, TI 
TIP4P, BAR 
T1P4P, FEP 


Experimental [28] 


-9.54 





"No long range van der Waals corrections 

TABLE 1: The hydration free energy of acetamide. (a) The Helmholtz hydration free energy AF for the current work, as 
calculated by the Bennett Acceptance Ratio, and the RMS hysteresis error. The A_F statistical errors are calculated by the 
bootstrap method, (b) Published values of the Gibbs free energy AG, obtained both computationally and experimentally. All 
computational results utilize the OPLS-AA force field for the solute acetamide. Also noted are the water model and free energy 
estimator (TI: Thermodynamic Integration; FEP: Free Energy Perturbation; BAR: Bennett Acceptance Ratio) 




FIG. 2: (a) The hysteresis error between neighboring replicas. Replica exchange effectively reduces the hysteresis error for 
replica pairs, (b) Block averages of the RMS hysteresis error, showing that the hysteresis error falls with increasing block 
size. Replica exchange increases the rate at which hysteresis error is lowered, thereby achieving the same magnitude error with 
simulations which are on average 4-8 times shorter. 



equilibrium ensemble, giving rise to a pronounced spike 
in C\. (Smaller shifts in C\ near A = 1 reflect electro- 
static effects and are not observed for simulations where 
Ac = 0, data not shown.) Thus, C\ profiles serve as use- 
ful probes for detecting large shifts in equilibrium ensem- 
bles. Regions where the equilibrium ensembles change 
most rapidly are the regions that contribute to the largest 
errors in free energy calculations. 



B. Optimal A Schedule for Free Energy 
Calculations 

For given computational resources, with the number of 
replicas and the simulation length fixed, the RMS hys- 



teresis error of a simulation may be decreased by opti- 
mizing the A schedule, or the distribution of A across the 
replicas. The swap probability gives the rate at which the 
average hysteresis error falls between two replicas, and in 
an optimized simulation it would be uniform across all 
replica pairs. In practice it is difficult to obtain the A 
schedule which makes the swap probability exactly uni- 
form, but reasonable approximations can be made by us- 
ing the linearized swap probability, given by Eq. (jlSp . 

First, it is necessary to perform some number of prelim- 
inary simulations to obtain C\ along a coarse A schedule. 
These initial simulations need not be as long as the final 
production runs, since C\ converges more quickly than 
6F and is more tolerant of error. With a rough estimate 
of Ca(A) in hand, the A schedule can be adjusted to en- 
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Mean Swap Probability 
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FIG. 3: The average swap probability between adjacent repli- 
cas and C\=v&r{dU/dX) evaluated for each replica (from the 
replica exchange simulation; simulation with no replica ex- 
change is not significantly different). Spikes in Cx indicate 
regions of low swap probability. 



Water Density in Vicinity of Acetamide 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



X value 

FIG. 4: Water density within 2.5A of the acetamide carbonyl 
carbon as A varies. The inset illustrates the position and size 
of the observation volume with respect to an ACE molecule. 
Density is normalized by the bulk density. As A increases, 
waters are expelled by the growing cavity. 



sure that the linear swap probability is uniform between 
all replicas. Alternatively, one might simply shift repli- 
cas from where C\ is small to where it is large. Both 
approaches are only approximate, and break down when 
the linear response assumption in Eq. (flSl) ceases to be 
valid. They may be applied iteratively as C\ is evaluated 
for new A schedules. 

The aim of an optimal A schedule is to place repli- 
cas close together in regions where the Cx profile shows 
spikes. This ensures reasonable swap probabilities and 
minimal hysteresis errors in regions that are problem- 



atic. Preliminary investigations show that even when the 
schedule is improved in an ad hoc manner, hysteresis as 
well as statistical errors decrease. 



C. Replica Exchange 

Replica exchange provides a Monte Carlo move which 
may allow a replica to access a distant part of its equi- 
librium ensemble in one step. It is no substitute for 
conformational exploration within a replica. This point, 
while obvious, must be emphasized in the context of the 
hysteresis error, which does not report on the quality of 
intra-replica sampling. As an extreme but illustrative 
case, consider a system of some number of frozen repli- 
cas, each with a different configuration, which undergo 
replica exchange moves but no conformational changes. 
With just a modest number of swaps, these configura- 
tions attain the probability distribution described by Eq. 
(fT3|) . and the hysteresis error is zero. The system has 
achieved inter-replica equilibrium, but the intra-replica 
probability distribution has not been obtained. In prac- 
tice, the majority of Monte Carlo moves must be within 
a replica. The optimal frequency of swap moves remains 
an open question, although preliminary simulations sug- 
gest that more frequent swaps reduce the hysteresis error 
more quickly. 



VI. SUMMARY AND CONCLUSION 

In a simulation of multiple replicas, each sampling the 
equilibrium ensemble of a different Hamiltonian, swap- 
ping configurations between replicas is a nonequilibrium 
work process. Accordingly, the work needed to perform 
such swaps has a distribution of values, as described by 
the Crooks fluctuation theorem. The hysteresis error en 
developed here measures how closely a given simulation 
reproduces these work distributions between a pair of 
replicas. 

The hysteresis error is particularly useful in the context 
of free energy calculations. It reports on the combined 
bias of the forward and reverse free energy perturbation 
techniques, and it measures how completely individual 
replicas sample their equilibrium ensemble. The RMS 
hysteresis error, which reports on en for the whole A 
schedule, may be decreased by running a longer simula- 
tion, employing replica exchange, utilizing an improved 
A schedule, or all of these approaches. 

The average swap probability is another useful measure 
and can be calculated whether or not replica exchange is 
employed. Since it determines the rate at which the hys- 
teresis error decreases with simulation length, the swap 
probability can be used to optimize the A schedule. With 
a uniform average swap probability the hysteresis error 
falls evenly between all replica pairs. This maximizes 
the efficiency of simulations with fixed computational re- 
sources, avoiding unnecessary replicas where the hystere- 
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sis is low and preventing excessive errors from regions 
where the hysteresis error is large. 

Furthermore, the swap probability, along with a re- 
lated measure C\ , yields insight into the microscopic be- 
havior of a system. The swap probability is low and C\ 
is large when the equilibrium ensemble changes rapidly 
with A. Slow convergence and bias errors in free energy 
calculations arise when there are spikes in the C\ profile 
along the A schedule, which coincides with large hystere- 
sis errors. 
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APPENDIX A: DERIVATIONS 

1. Fluctuation theorem derivation 

We derive the Crooks fluctuation theorem ^ in the 
context of instantaneously switching Aq — > Ai (forward) 
and Ai — » Aq (reverse). Expanding the ratio po/pi with 
([3]) for an arbitrary configuration F, 



po(r) 
Pi(r) 



and similarly. 



exp [p{Fo-F^)-P{Uo-U{)] 



Pi(r) 

po(r) 



= exp[/3VF«(r)] 



(Ala) 



(Alb) 



where the definitions of work @ and dissipated work ([7]) 
were used. 

We integrate pi from (jAlap over all configurations, but 
consider contributions only from those F for which the 
forward dissipated work work value takes on a specific 
value, Wd' 



I 



drpo{r)eM-f3wE{TmpWD - pwE{r)] 

dTpi{T)S[f3WD - PWEir)]. (A2) 



Since, from (jAlaP and (lAlbp . 

wEir) = -wg{T) 



21) becomes, 

dTpo{T)eM-PWE{T)]6[pWD ~ PWE{T)] 

drpi{T)S[pWD + PWEi^)]- (A3) 



We define P^{Wd) as the probability of observing a 
given dissipated work value in the forward switching pro- 
cess, and it can be expressed as an integral over all con- 
figurations which yield this value, 

P^(Wd) - J drpoiT)6[pWD - PWEiT)] (A4a) 

Likewise, the probability of observing a given disspated 
work value in the reverse switching process is, 

P"{Wd) = / drp,{r)6[(3WD - (3WE{T)] (A4b) 



With these definitions, (jA3p may be written as, 

exp{-f3WD)P^{PWD) = P^i-pWo), 
which is equivalent to ([T|). 

2. Fluctuation theorem and hysteresis error 

The relationship between some arbitrary deviation of a 
simulation from the Crooks fluctuation theorem and the 
hysteresis error is derived by first rewriting Eq. ([8]) as, 

P*ni-f3WD)eMf3WD) = P*Al3WD)eM-(3e*FT)- (A5) 

Inserting the SF^^p definition (|5bp into the definition of 
the hysteresis error ^ , expanding the reverse work with 
(j7bp and using the SFEep estimate for 5F, 

en - SFEep - P'^HeM-PW^)Yi, 

= SFEep - In [{eM-PWE))l eMP^FEEp)] , 

We now expand the estimated ensemble average as an 
integral over all values of f3WE, with P^ the normalized 
histogram of (3WE obtained from a simulation. 



~ (3-^11 



d[pwg]PAPwE)eM-Pwg) 



As pWl^ is a dummy variable, we change it to —fSWo 
en = In ' 



d[l3WD]PR{-l3WD) eMPWo) 



where we implicitly multiplied the integrand by —1 to 
preserve the limits of integration. With (|A5p the above 
can be written as, 



EH = -/3-iln 



-f-oo 



d{l3WD)PF{PWD)exp{-f3 



which reduces to 
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3. Inter-replica equilibrium and hysteresis error 

We can relate an small arbitrary error in the calcu- 
lated distribution to the hysteresis error by consid- 
ering a small error /9j(ro,ri) in the otherwise correctly 
estimated p'^. Rewriting (fT3| . 



we integrate over all configuration pairs and rewrite 

AUsroap with ((nil, 



dTodTip'j^ + J dTodTip^ 

= J drodripo(ro)pi(ri)exp[-/3W^^(ro)] 

X y"dripi(ri)exp[-/3M/«(ri)]. (A6) 



With the sampling error contained in pe, the p^ term 
(expanded with (|12bp ) is identically one. Taking the log- 
arithm and dividing by /3, (|A6|) becomes, 



p-^ In 



1 



dTodTip^ 



^SF^EP-^F^EP. (A7) 



where we have used the SFfep definitions ([5]). With the 
approximation ln(l+a::) ~ x for small x and the definition 
of Cif dS]), we obtain Eq. ([M]) . 



4. Linearized average swap probability 

Here we consider the average Fermi swap probability 
between two replicas whose A parameters differ by a small 
amount, 5 (written as 5\ in the text). For convenience 
we define 

P = pALfg^,fipj 

= f3[UsiTo) - UoiTo) + Uo{Ts) - UsiTg)], 

where Fq and F^ are configurations drawn from the equi- 
librium distributions Uq and Us parameterized by Aq and 
Aq + S, respectively. We expand Us as a Taylor series 
about Aq, 



Us{T) = [/o(F) + 5Vo{T) + -WoiT) + 0(6^), 



with 



V„ 



Wo 



dU_ 

d^U 



5A2 



A=Ao 



p can then be written as, 



p^pS[VoiTo)-VoiTs)] 



[Wo{To)-WoiTs)]. 



Note that p is small {0{5)); thus, with the identities, 

exp(a::) = 1 + x + x'^ /2 + (A8a) 



1 + x 



= 1- x + X'' 



(A8b) 



we may write the Fermi swap probability between con- 
figurations Fq and F^ as, 



1 



Pswap — 



1 + cxp p ' 

1 ^ 1 



2 Vi + m/2 + mV4 + 0(m3)/ ' 
^ ^1 - (p/2 + pV4) + (p/2 + pV4)2 + 0(p3)] , 



2 

1 1 

2 ~ 4 



p + 0(p3). 



The average swap probability is the ensemble average 
over all configuration pairs, 

{{Pswap)Q)s = ^ - ^((m)o)5, 



-PS{Vo)s~^{Wo)s 



(A9) 



To evaluate {■)s, we first obtain Qs, the partition function 
at (Xo + S): 

Qs = J dTcxpi^pUs), 

= I dTcxp{-l3Uo) [1 - mo + 0{d^)] , 
= Qo[l- f3d{Vo)o + 0{5^)] , 
and its reciprocal, 

Qj' = Qo' [l + l3S{Vo)o + 0{S^)] . 

We can now evaluate {Vq)s and {Wo)s, retaining only 
terms which will remain 0{d^) or larger in (jA9p : 



{Vo)s = Qj' J dVeM-PUsWo, 

= Qo' {I + f36 {Vo)o) I dF(l-/3(51/o)exp(-/3C/o)yo, 



= (l + /3(5(Vb)o)((K))o-/?'5(K)'>o), 

and 

{Wo)s = Qj' J dTeM-PUs)Wo, 

= Qo\l + 0{S)) J dFexp(-/?f/o)Wo[l-0(<S)], 
= {Wo)o + OiS). 
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Finally, (IA9p becomes, 
1 (3^6^ 

( {Pswap )o)s 



{{V,\-{Vo)l)+0{5'), (AlO) 



2 4 

equivalent to Eq. (HH]), which is valid for small 6. 



APPENDIX B: Ulj AND Uc FUNCTIONAL 
FORMS 

The functional forms of both the Coulomb and 
Lennard-Jones potentials were developed for this work 
based on three criteria: 

1. Configurations where the solute and solvent overlap 
may be observed for A = 0. For such configurations, 
we require: 

• That swaps be permitted with reasonable fre- 
quency for small A (e.g. A — 0.1). 

• That swap probabilities falls off quickly there- 
after; in particular, we wish to avoid the situ- 
ation where the swap probability declines only 
very near A = 1.0. 

2. We require that dU /dX is not always zero for A = 
to avoid complications with the Thermodynamic 
Integration (TI) estimator. While, we do not report 
results using TI in this work, we wish to construct 
a A schedule that works with all estimators. 

3. In this work, A^j — Ac- Therefore, Lennard-Jones 
repulsion must dominate Coulombic attraction at 
very small atomic separations. 

While various ways to scale the potential have been dis- 
cussed in the literature [1^, [s^, HH , none of these satisfied 



all of our requirements. It should be noted that condi- 
tion 3 is somewhat arbitrary, and more common scaled 
potentials may be used if the insertion process scales the 
Lennard-Jones prior to the Coulomb potential. 

a. Coulomb scaling We employ a modified version of 
the linear soft-core scaling [13] ; for two atoms of charges 
Qi and Qj distance r apart, the potential energy is Ac as. 



Uc{r, Ac) = Ac 



ac(l - Ac) + r 



(Bl) 



ac controls the "soft core" term, and for small Ac im- 
poses a minimum effective atomic separation, ac — 1.5 A 
for all simulations in this work. 

b. Lennard-Jones scaling The Lennard-Jones po- 
tential between two particles may be written generally 



ULj{r,XLj)^BAiA-l) 
where, for unsealed Lennard-Jones, 



(B2) 



^(r) = (^) , B = Ae. 

Simple linear scaling by Xlj of the Lennard-Jones po- 
tential is known to be unsatisfactory, and a number of 
alternate forms have been introduced. We developed the 
exponential soft-core. 



A{r,XLj) - 1/ 
B{Xlj) = 4e. 



aLjil-XLjf 



(- 



1 



, (B3a) 
(B3b) 



with a = 4, fc = 1 and a^j ~ O.bA. The precise position 
along the A coordinate of the swap probability trough 
(see Fig. [3]) is specific to this Lennard-Jones potential. 
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